% code to compute MT trfn from raw ts files !!!
% written for svs
% basically to yeild tfs at all FC's
% clears unneceesary vriables after
% invoking readrts reads the header
% compute cal. table according to gain settings etc;
% computes cross and auto spectra & average it 
% MT app rest & phse are now computed
% - Bug report: fixed the conjugate problem
% latest date 6.10.5

function[] = eqspec(),
readrts;

Notchi = NotchF(2:5);		% NOTE ! there may be problems 
IsNotch = Filt(2:5);		%        if notch do not have right dimension..
hpb1 = Filt(1);			% NOTE ! high pass filter is always disabled
prob = [Exd Eyd Eyd Eyd Eyd];	% 	 in calib05.m
gain = [Exa Eya Hxa Hya Hza];
clear Bden Record3 Bdend Bddm Count Record4 hpf Bdst Record5;
clear Bdstart MaxRecord ans nt Filt NBlocks array nlines Events ;
clear NotchF dummy Record1 factor Channel;
BandFrq = [8192 256 8 1/4 1/128 1/4096]; % Here the upper and lower limit of frequencies of each band is given
RadSize = [0.125 0.25 0.5 0.75];
choiseR = 1;
block = 1024;
TLFreq = logspace(log10(BandFrq(nband)),log10(BandFrq(nband+1)),20);
TLFreq([1 end]) = [];
notch = Notchi(find(IsNotch>0)); 
PAI = 3.1415926535897;

ProcDef = setdefault(nband);

[TLine,TLRad,TLFreq,TBand] =  SetTLine(block,RadSize(choiseR),TLFreq);
ProcDef.TLFreq1 = TLFreq;
for i = 1:nchannel,
[cplx,frq1] = calib05(i,nband,block,prob(i),notch,hpb1,gain(i));
CalData(i,:) = cplx;
end;

%------------native declarations---------------

WI = hanning(block);
len = length(matrix)/block;
start = 0; %change len and start to obtain stacked spectra of desired interval
len = 60;
st = start*block + 1;
en = (start+1)*block;
nchannel = 5;


nfrq = length(frq1);
HxHx = zeros([nfrq,1]);
HyHy = HxHx;
HxHy = HxHx;
HyHx = HxHx;
ExHx = HxHx;
EyHx = HxHx;
ExHy = HxHx;
EyHy = HxHx;
ExEx = HxHx;
EyEy = HxHx;

icount = 0;

for i = (start+1):len,
icount = icount+1;
%fft & calib
for u = 1:5,
	K = CorTre(matrix(st:en,u),0);
	Channel1(:,u) = fft(K.*WI)*(2/block);
	Channel(:,u) = Channel1(2:(block/2)+1,u).*conj(CalData(u,1:block/2)');
end; 	
for ii = 1:length(TLFreq),
   par = parzen1(TLRad(ii)); % Bit inefficient (as it is called many times to do the same job)
   SPMatrix(ii,i,:,:) = CompuSPM1(Channel,par,TLine(ii),TLRad(ii));
end;

% auto and cross spectra---------------------------
   
   HxHx   = HxHx + Channel(:,3).*conj(Channel(:,3));
   HyHy   = HyHy + Channel(:,4).*conj(Channel(:,4));
   HxHy   = HxHy + Channel(:,3).*conj(Channel(:,4));
   HyHx   = HyHx + Channel(:,4).*conj(Channel(:,3));
   ExHx   = ExHx + Channel(:,1).*conj(Channel(:,3));
   EyHx   = EyHx + Channel(:,2).*conj(Channel(:,3));
   ExHy   = ExHy + Channel(:,1).*conj(Channel(:,4));
   EyHy   = EyHy + Channel(:,2).*conj(Channel(:,4));
   ExEx   = ExEx + Channel(:,1).*conj(Channel(:,1));
   EyEy   = EyEy + Channel(:,2).*conj(Channel(:,2));

%------------------------------------------------

%clear Channel;

st = en + 1;
en = (i+1)*block;
end;


HyHy =HyHy/icount;
HxHy =HxHy/icount;
HyHx =HyHx/icount;
ExHx =ExHx/icount;
EyHx =EyHx/icount;
ExHy =ExHy/icount;
EyHy =EyHy/icount;
HxHx =HxHx/icount;
ExEx = ExEx/icount;
EyEy = EyEy/icount;

% 4.1.06 Whne compared with Mappros Cal Stack spectra (1024)
% a factor to of (2/100) became necessary to mak ethe amplitudes similar
% (no idea about the orgin of the factor)
% So ExE = ExEx*(2/100), HyHy = HyHy*(2/100) etc


Zxx = HyHy.*ExHx - HyHx.*ExHy;
Zyy = HxHx.*EyHy - HxHy.*EyHx;

Zxy = HxHx.*ExHy - HxHy.*ExHx;
Zyx = HyHy.*EyHx - HyHx.*EyHy;
Den = HxHx.*HyHy - HxHy.*HyHx;



for i = 1:nfrq,
   if abs(Den(i)>0)&abs(Zxy(i))>0&abs(Zyx(i))>0,
      ZXY(i)=Zxy(i)/Den(i);
      ZYX(i)=Zyx(i)/Den(i);
      ZXX(i)=Zxx(i)/Den(i);
      ZYY(i)=Zyy(i)/Den(i);
   end;
end;




for i = 1:nfrq,
   Phxy1(i) = atan(imag(ZXY(i))/real(ZXY(i))) * 180 / 3.1415926535897;
   Rhxy1(i) = 0.2 * (1/frq1(i)) * (ZXY(i) * conj(ZXY(i)));
   Phyx1(i) = atan(imag(ZYX(i))/real(ZYX(i))) * 180 / 3.1415926535897;
   Rhyx1(i) = 0.2 * (1/frq1(i)) * (ZYX(i) * conj(ZYX(i)));
end;


%--------------------------------------------------------------
 T = 1./frq1;
 
subplot(2,1,1);
loglog(T,Rhxy1,'r.',T,Rhyx1,'b.');
ylabel('ohm-m');
xlabel('period (s)');
grid;
subplot(2,1,2);
semilogx(T,Phxy1,'k.',T,Phyx1,'k.');
xlabel('period (s)');
ylabel('degrees');
title(file);
grid;

figure(2);
SPM2  = AverageSPM(SPMatrix,ProcDef);
[tf2]=tf_dn(SPM2);
[Rhxy Rhyx Phxy Phyx] = dispresph(tf2,ProcDef);




% computes target frequency, target lines and radius 
% for current block size;
% note Parzen window sizes are hard written here
% part of Mt processing codes
% usage [TLine,TLRad,TLFreq,TBand] =  SetTLine(block)
% latest date 07.11.2001
% 

function[TLine,TLRad,TLFreq,TBand] =  SetTLine(block,RadSize,TLFreq),

	DeltaFreq = [32.0 1.0 1.0/32.0 1.0/1024.0 (1.0/32768.0)/0.25];
 	MinFreq      = 1.0 / 4097.0;
	MaxFreq      = 8193.0;



for Index = 1:length(TLFreq),

    if TLFreq(Index) < 8.0/1024.0,
      TLFreq(Index) =  floor(TLFreq(Index) / DeltaFreq(5) + 0.5)*DeltaFreq(5);
	TBand(Index) = 5;
    elseif TLFreq(Index) < (8.0/32.0),
      TLFreq(Index) =  floor(TLFreq(Index) / DeltaFreq(4) + 0.5)*DeltaFreq(4);
	TBand(Index) = 4;
    elseif TLFreq(Index)   < 8.0 ,
      TLFreq(Index) =  floor(TLFreq(Index) / DeltaFreq(3) + 0.5)*DeltaFreq(3);
	TBand(Index) = 3;
    elseif TLFreq(Index) < 256.0,
      TLFreq(Index) =  floor(TLFreq(Index) / DeltaFreq(2) + 0.5)*DeltaFreq(2);
	TBand(Index) = 2;
	else,
      TLFreq(Index) =  floor(TLFreq(Index) / DeltaFreq(1) + 0.5)*DeltaFreq(1);
	TBand(Index) = 1;
    end;


    if TLFreq(Index) <= MinFreq,
      TLFreq(Index)  = MinFreq;
    elseif TLFreq(Index) >= MaxFreq,
      TLFreq(Index)  = MaxFreq;
    end;


%    -  set targetlines and band                            
%    -  targetlines are set for block size of 1024 points 

    if TLFreq(Index) < 8.0/1024.0,
      TLine(Index) =  floor(TLFreq(Index) / DeltaFreq(5) + 0.5);
    elseif TLFreq(Index) < (8.0/32.0),
      TLine(Index) =  floor(TLFreq(Index) / DeltaFreq(4) + 0.5);
    elseif TLFreq(Index)   < 8.0 ,
      TLine(Index) =  floor(TLFreq(Index) / DeltaFreq(3) + 0.5);
    elseif TLFreq(Index) < 256.0,
      TLine(Index) =  floor(TLFreq(Index) / DeltaFreq(2) + 0.5);
	else,
      TLine(Index) =  floor(TLFreq(Index) / DeltaFreq(1) + 0.5);
    end;


%   -- set correct targetline in respect to block 


	if block == 128,
	TLine(Index) = round((TLine(Index) +1 ) / 8);
	elseif block == 256,
	TLine(Index) = round((TLine(Index) +1) / 4);
	elseif block == 512,
	TLine(Index) = round((TLine(Index)) / 2);
	elseif block == 1024,
	TLine(Index) = TLine(Index);
	elseif block == 2048,
	TLine(Index) = round((TLine(Index) +1) * 2);
	elseif block == 4096,
	TLine(Index) = round((TLine(Index) +1) * 4);
end;



%------  set radius  ---------------------------------------


TLRad(Index) = floor(TLine(Index) * RadSize + 0.5 );

%----- check radius boundaries 

 %   if TLRad(Index) < MinParzenRad,
 %     TLRad(Index) = MinParzenRad;
%    elseif TLRad(Index) > MaxParzenRad,
%      TLRad(Index) = MaxParzenRad;
%    end;

end;


function[SPM] = CompuSPM1(Channel, par,TL,TR),

a = size(Channel);
nchannel = a(2);

for i = 1:nchannel,
   for j = 1:i,
talkal = sum(Channel(TL-(TR-1):TL+(TR-1),j).*conj(Channel(TL-(TR-1):TL+(TR-1),i)).*par); % summing w.r.t parzen window%
      if i ~= j,
       SPM(i,j) = real(talkal);
       SPM(j,i) = (-1)*imag(talkal); 
      elseif i == j,
       SPM(i,j) = talkal;
     end;
  end;
end;

function[PF] = parzen1(radius),
    if radius > 1,

%---- compute norm-factor for parzening 

      nFactor = 1.0; 
		PF(1) = 1.0;
      for line = 2:radius,

        u = (line/(radius+1))*pi;
        v = (sin(u)/u)^2.0;
        PF(line) = v;
        nFactor = nFactor+(2*v*v);
      end;
      nFactor =sqrt(nFactor);

%      -- parzening line 1 

      PF(1) = 1.0/nFactor;

%     -- parzening line 2 to radius      
	for line = 2:radius,
        PF(line) = PF(line)/nFactor;
     end;
     ParInt = sum(PF(2:radius))*2 + PF(1);
PF = PF/ParInt;
else,
   PF = 1;
   end;
   PF= [PF(radius:-1:2) PF]';
   
   
   %----trying to make parzen average a robust est 25.5.3''


